Hydrology in R demonstration airGR - EGU 24/04/2017

Session lead: Katie Smith (Centre for Ecology & Hydrology) k.a.smith@ceh.ac.uk

This is an demonstration of how to use the airGR package of hydrological models in R, as well as how to plot interactive timeseries graphs with the dygraphs package.

First we need to install some packages

library(xts)

package <U+393C><U+3E31>xts<U+393C><U+3E32> was built under R version 3.3.3Loading required package: zoo

Attaching package: <U+393C><U+3E31>zoo<U+393C><U+3E32>

The following objects are masked from <U+393C><U+3E31>package:base<U+393C><U+3E32>:

    as.Date, as.Date.numeric

library(dygraphs)

package <U+393C><U+3E31>dygraphs<U+393C><U+3E32> was built under R version 3.3.3

Now we’ll load in some observational flow data from the River Thames (naturalised) in England - with thanks to the National River Flow Archive: http://nrfa.ceh.ac.uk/data/search

observed_data <- read.csv("~/Easter_Work/GR_teaching/Qobs_390010.csv")
observed_data

observed_data$DATE <- strptime(observed_data$DATE, format = "%d/%m/%Y")

in order to plot this as an interactive dygraph we need to change it to xts format

observed_data_xts <- as.xts(observed_data$Qobs, order.by = observed_data$DATE)
dygraph(observed_data_xts)%>%
  dyOptions()%>%
  dyRangeSelector()

Now lets read in some precipitation data - this is from CEH-GEAR: https://data.gov.uk/dataset/gridded-estimates-of-daily-and-monthly-areal-rainfall-for-the-united-kingdom-1890-2012-ceh-gear

precip_data <- read.csv("~/Easter_Work/GR_Teaching/rain_1961_2014_390010.csv")
precip_data

and some potential evapotranspiration data - this is from CHESS-PE:https://data.gov.uk/dataset/climate-hydrology-and-ecology-research-support-system-potential-evapotranspiration-dataset-for-1

PET_data <- read.csv("~/Easter_Work/GR_Teaching/CHESS_PET_1961_2015_390010.csv")
PET_data

Note that our starting dates are not the same as our observational data, so we need to make a dataframe that matches the dates up. There are a lot of ways to do this. First we’ll convert them to the same date format.

precip_data$DATE <- strptime(precip_data$DATE, "%Y-%m-%d")
PET_data$DATE <- strptime(PET_data$DATE, "%Y-%m-%d")

now we’ll find the common period

first_date <- max(observed_data$DATE[1], precip_data$DATE[1], PET_data$DATE[1])
last_date <- min(observed_data$DATE[length(observed_data$DATE)], precip_data$DATE[length(precip_data$DATE)], PET_data$DATE[length(PET_data$DATE)])

and make a data frame of that length

# make an empty data frame
thames_data <- as.data.frame(matrix(NA,nrow=as.numeric((last_date-first_date)+1), ncol=4))
colnames(thames_data) <-c ("date","PET","precip","obs")
# make the date timeseries
thames_data$date <- seq(first_date, last_date, by="days")
# populate the data frame with the data
thames_data$obs <- observed_data$Qobs[which(observed_data$DATE==thames_data$date[1]):which(observed_data$DATE==thames_data$date[length(thames_data$date)])]
thames_data$precip <- precip_data$Mean_rainfall[which(precip_data$DATE==thames_data$date[1]):which(precip_data$DATE==thames_data$date[length(thames_data$date)])]
thames_data$PET <- PET_data$PET[which(PET_data$DATE==thames_data$date[1]):which(PET_data$DATE==thames_data$date[length(thames_data$date)])]

plot the observed streamflow with the precipitation data

#  convert the observed discharge to runoff (so its in the same units as the precip)
# divide by catchment area (m2) and mulitply by 86.4
thames_data$obs <- (thames_data$obs/9948.0)*86.4
thames_data_xts <- as.xts(thames_data[,3:4], order.by=thames_data$date)
# initiate the dygraph
dygraph(thames_data_xts, main = "Naturalised Runoff and Precipitation Observations for the Thames at Kingston")%>%
# define the first axis  
dyAxis(dygraph = graphOut, name = "y", label = "runoff (mm/day)",
       valueRange = range(thames_data_xts[, "obs"],
                          na.rm = TRUE)* c(0.01, 1.59))%>%
# define the second axis
dyAxis(dygraph = graphOut, name = "y2", label = "precip (mm/day)",
                   valueRange = rev(range(thames_data_xts[, "precip"], 
                   na.rm = TRUE)* c(0.01, 2.99)))%>%
# plot the data
dySeries("obs",axis = 'y')%>%
dySeries("precip", axis = 'y2', stepPlot = TRUE,
         fillGraph = TRUE)%>%
dyOptions(colors = RColorBrewer::brewer.pal(3,"Set1")[3:1]) %>%
dyRangeSelector()

OK, enough messing with data, lets do some modelling

see this website for a good guide through the model: https://odelaigue.github.io/airGR/tutorial_1_getting_started.html

import the GR package

require(airGR, quietly=TRUE)

package <U+393C><U+3E31>airGR<U+393C><U+3E32> was built under R version 3.3.3

prepare the input data in the correct format

BasinObs <- thames_data
colnames(BasinObs) <- c('DatesR','E','P', 'Qobs')

create the InputsModel object - this defines which model we want to run, and defines the variables for the models input data

InputsModel <- CreateInputsModel(FUN_MOD = RunModel_GR4J,DatesR = BasinObs$DatesR,
                                   Precip = BasinObs$P,PotEvap = BasinObs$E)
str(InputsModel)

List of 3
 $ DatesR : POSIXlt[1:19723], format: "1961-01-01 00:00:00" "1961-01-02 00:00:00" ...
 $ Precip : num [1:19723] 8.6376 4.571 2.1015 0.0201 7.1859 ...
 $ PotEvap: num [1:19723] 0.367 0.353 0.412 0.32 0.612 ...
 - attr(*, "class")= chr [1:3] "InputsModel" "daily" "GR"

# note NA values of precip and PET are NOT ALLOWED

create the RunOptions object - this defines options for the RunModel_GR4J function

# first define the period to run the model over
Ind_Run <- seq(which(BasinObs$DatesR=="1981-01-01"),
             which(BasinObs$DatesR=="2014-12-31"))
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
                               InputsModel = InputsModel,
                               IndPeriod_Run = Ind_Run,
                               IndPeriod_WarmUp = NULL)

     Model warm-up period not defined -> default configuration used 
        The year preceding the run period is used 
     Model states initialisation not defined -> default configuration used 

str(RunOptions)

List of 6
 $ IndPeriod_WarmUp: int [1:365] 6941 6942 6943 6944 6945 6946 6947 6948 6949 6950 ...
 $ IndPeriod_Run   : int [1:12418] 7306 7307 7308 7309 7310 7311 7312 7313 7314 7315 ...
 $ IniStates       : num [1:67] 0 0 0 0 0 0 0 0 0 0 ...
 $ IniResLevels    : num [1:2] 0.3 0.5
 $ Outputs_Cal     : chr "Qsim"
 $ Outputs_Sim     : chr [1:16] "DatesR" "PotEvap" "Precip" "Prod" ...
 - attr(*, "class")= chr [1:3] "RunOptions" "GR" "daily"

create the InputsCrit object - define the error metric (choose from RMSE, NSE, KGE or modified KGE (KGE2))

InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_NSE,
                               InputsModel = InputsModel,
                               RunOptions = RunOptions,
                               Qobs = BasinObs$Qobs[Ind_Run])
str(InputsCrit)

List of 5
 $ BoolCrit  : logi [1:12418] TRUE TRUE TRUE TRUE TRUE TRUE ...
 $ Qobs      : num [1:12418] 0.635 0.623 0.587 0.571 0.566 ...
 $ transfo   : chr ""
 $ Ind_zeroes: NULL
 $ epsilon   : NULL
 - attr(*, "class")= chr "InputsCrit"

create the CalibOptions object - choose the calibration algorithm

CalibOptions <- CreateCalibOptions(FUN_MOD = RunModel_GR4J,
                                   FUN_CALIB = Calibration_Michel)
str(CalibOptions)

List of 3
 $ FixedParam       : logi [1:4] NA NA NA NA
 $ SearchRanges     : num [1:2, 1:4] 4.59e-05 2.18e+04 -1.09e+04 1.09e+04 4.59e-05 ...
 $ StartParamDistrib: num [1:3, 1:4] 169.017 247.151 432.681 -2.376 -0.649 ...
 - attr(*, "class")= chr [1:3] "CalibOptions" "GR4J" "HBAN"

run the calibration

OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, 
                                   RunOptions = RunOptions,
                                   InputsCrit = InputsCrit,
                                   CalibOptions = CalibOptions,
                                   FUN_MOD = RunModel_GR4J,
                                   FUN_CRIT = ErrorCrit_NSE)

Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
     Screening completed (81 runs)
         Param =  432.681 ,   -0.649 ,   83.096 ,    2.384
         Crit NSE[Q]       = 0.8923
Steepest-descent local search in progress
     Calibration completed (19 iterations, 216 runs)
         Param =  607.894 ,   -0.734 ,   87.357 ,    2.315
         Crit NSE[Q]       = 0.9246

NSE of 0.9246 - not bad at all!

define the parameters found by the calibration routine

Param <- OutputsCalib$ParamFinalR
Param

[1] 607.8936811  -0.7336304  87.3567230   2.3153153

RUN THE MODEL!

OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
                             RunOptions = RunOptions,
                             Param= Param)
str(OutputsModel)

List of 16
 $ DatesR  : POSIXlt[1:12418], format: "1981-01-01 00:00:00" "1981-01-02 00:00:00" ...
 $ PotEvap : num [1:12418] 1.347 0.38 0.795 0.57 0.454 ...
 $ Precip  : num [1:12418] 0.0972 0.1211 0.0523 0.1147 0.3248 ...
 $ Prod    : num [1:12418] 352 351 350 350 349 ...
 $ AE      : num [1:12418] 1.127 0.334 0.663 0.488 0.43 ...
 $ Perc    : num [1:12418] 0.387 0.383 0.378 0.374 0.372 ...
 $ PR      : num [1:12418] 0.387 0.383 0.378 0.374 0.372 ...
 $ Q9      : num [1:12418] 0.355 0.35 0.345 0.341 0.338 ...
 $ Q1      : num [1:12418] 0.0398 0.0393 0.0387 0.0382 0.0378 ...
 $ Rout    : num [1:12418] 42.2 42 41.7 41.4 41.2 ...
 $ Exch    : num [1:12418] -0.0592 -0.0577 -0.0564 -0.0551 -0.0539 ...
 $ AExch   : num [1:12418] -0.099 -0.097 -0.0951 -0.0933 -0.0917 ...
 $ QR      : num [1:12418] 0.599 0.578 0.559 0.542 0.526 ...
 $ QD      : num [1:12418] 0 0 0 0 0 ...
 $ Qsim    : num [1:12418] 0.599 0.578 0.559 0.542 0.526 ...
 $ StateEnd: num [1:67] 374.4 45.5 NA NA NA ...
 - attr(*, "class")= chr [1:3] "OutputsModel" "daily" "GR"

use the inbuilt plot function to look at the results

plot(OutputsModel, Qobs=BasinObs$Qobs[Ind_Run])

looking good - but we’ve got some discrepancy at the low flows end. NSE is notorious for this, it is based on the square of the flows, so over-weights the calibration to the high flows. I wonder if the modified KGE can do any better?

# make a few changes to the calibration criteria
InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2,
                               InputsModel = InputsModel,
                               RunOptions = RunOptions,
                               Qobs = BasinObs$Qobs[Ind_Run])
# rerun the calibration
OutputsCalib <- Calibration_Michel(InputsModel = InputsModel, 
                                   RunOptions = RunOptions,
                                   InputsCrit = InputsCrit,
                                   CalibOptions = CalibOptions,
                                   FUN_MOD = RunModel_GR4J,
                                   FUN_CRIT = ErrorCrit_KGE2)

Grid-Screening in progress (0% 20% 40% 60% 80% 100%)
     Screening completed (81 runs)
         Param =  432.681 ,   -0.649 ,   83.096 ,    2.384
         Crit KGE'[Q]      = 0.8589
Steepest-descent local search in progress
     Calibration completed (51 iterations, 502 runs)
         Param =  598.748 ,   -0.690 ,   93.679 ,    2.297
         Crit KGE'[Q]      = 0.9621

# redefine the parameters
Param <- OutputsCalib$ParamFinalR
# rerun the model
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
                             RunOptions = RunOptions,
                             Param= Param)
# plot again
plot(OutputsModel, Qobs=BasinObs$Qobs[Ind_Run])

not much different. Oh well, we can be happy with either of those metric scores. - pause for thought - which parameter set would you choose to use?!

Let’s do some validation

# go back to the beginning, redefine the period to run on (the period we haven't used for calibration, minus the first year needed for warm up)
Ind_Run <- seq(which(BasinObs$DatesR=="1962-01-01"),
             which(BasinObs$DatesR=="1980-12-31"))
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
                               InputsModel = InputsModel,
                               IndPeriod_Run = Ind_Run,
                               IndPeriod_WarmUp = NULL)

     Model warm-up period not defined -> default configuration used 
        The year preceding the run period is used 
     Model states initialisation not defined -> default configuration used 

InputsCrit <- CreateInputsCrit(FUN_CRIT = ErrorCrit_KGE2,
                               InputsModel = InputsModel,
                               RunOptions = RunOptions,
                               Qobs = BasinObs$Qobs[Ind_Run])
Param <- OutputsCalib$ParamFinalR
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
                             RunOptions = RunOptions,
                             Param= Param)
OutputsCrit <- ErrorCrit_KGE2(InputsCrit = InputsCrit, 
                              OutputsModel = OutputsModel)

Crit. KGE'[Q] = 0.9039
    SubCrit. KGE'[Q] cor(sim, obs, "pearson") = 0.9466 
    SubCrit. KGE'[Q] sd(sim)/sd(obs)          = 0.9253 
    SubCrit. KGE'[Q] mean(sim)/mean(obs)      = 1.0282 

slightly worse than the calibration period (0.9621) but not bad at all

finally, lets run the model for the whole time period and plot a dygraph so we can look at the timeseries in more detail

Ind_Run <- seq(which(BasinObs$DatesR=="1962-01-01"),
             which(BasinObs$DatesR=="2014-12-31"))
RunOptions <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
                               InputsModel = InputsModel,
                               IndPeriod_Run = Ind_Run,
                               IndPeriod_WarmUp = NULL)

     Model warm-up period not defined -> default configuration used 
        The year preceding the run period is used 
     Model states initialisation not defined -> default configuration used 

Param <- OutputsCalib$ParamFinalR
OutputsModel <- RunModel_GR4J(InputsModel = InputsModel,
                             RunOptions = RunOptions,
                             Param= Param)
plot_output_data <- as.data.frame(matrix(NA, ncol = 3,
                                         nrow = length(OutputsModel$DatesR)))
colnames(plot_output_data) <- c("Date", "Qsim", "Qobs")
plot_output_data$Date <- OutputsModel$DatesR
plot_output_data$Qsim <- OutputsModel$Qsim
plot_output_data$Qobs <- BasinObs$Qobs[Ind_Run]
plot_output_data_xts <- as.xts(plot_output_data, order.by = plot_output_data$Date)
library(RColorBrewer)
dygraph(plot_output_data_xts, main="Observed and Simulated Runoff for the Thames at Kingston (Naturalised)")%>%
  dyOptions(colors = RColorBrewer::brewer.pal(3,"Set1")[3:1])%>%
  dyAxis("y", label="Runoff (mm/day)")%>%
  dyRangeSelector()